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' We perform a numerical free evolution of a selfgravitating, spherically symmetric scalar 

pL^ ' field satisfying the wave equation. The evolution equations can be written in a very simple 

form and are symmetric hyperbolic in Eddington-Finkelstein coordinates. The simplicity 
of the system allow to display and deal with the typical gauge instability present in these 
coordinates. The numerical evolution is performed with a standard method of lines fourth 
^ . order in space and time. The time algorithm is Runge-Kutta while the space discrete deriva- 

' five is symmetric (non-dissipative) . The constraints are preserved under evolution (within 

, numerical errors) and we are able to reproduce several known results. 

o 

O ; 1 Introduction 



In the past decade, due to the results of several numerical experiments, there has been a consid- 
O ■ erable understanding of the phenomena of the collapse of a spherically symmetric, selfgravitating 

scalar field. In particular universal scaling properties of the final mass of the black hole were 



^ ■ discovered and then explained. For a review and references see The present work is not aimed 

to continue that line, but rather to study instability problems which are recurrent in many nu- 
merical evolutions near black holes when those evolutions are completely free, that is when the 
k> \ constraints are not solved for, but remain being satisfied just because the evolution equations are 

''^ ■ solved for. There are few cases in which free evolution has been successful, some of them, H, ^, ^ 

I use light cone coordinates and in fact some of the equations being solved could be considered 

as constraints. Others, ^, uses a symmetric hyperbolic formulation and use the freedom still 
available in that setting to suppress the main instability found there and so reach a stable prop- 
agation, alas first order numerical dissipation is used in that scheme. Finally there are conformal 
space-time methods ^ where no instability seems to be present, nor should be very relevant, 
for only short time evolutions suffices to cover the whole future domain of dependence. 

In this work we use Eddington-Finkelstein coordinates and are able to choose equations and 
variables so that the final system of equations is remarkably simple and manifestly symmetric 
hyperbolic. Such simplicity allows as to study the nature of the instability and proceed to avoid 
it. 

In the next section we introduce the equations and discuss which evolution system, among 
all possible equivalent ones, is suitable for free evolution. It turns out that two out of the five 
characteristics of the system can be specified at will by modifying the evolution equations. This 
is done by adding terms proportional to the constraints. These characteristics can even become 

*Work supported by Antorchas, CONICOR, CONICET and Se.CyT, UNC 
tPartially supported by AIT, Cordoba, Argentina. 
tMember of CONICET. 
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complex, signaling the failure of the evolution equations to be hyperbolic, and so to have a 
well posed initial value formulation. The change on the characteristics of the evolution system 
reflects on an identical change on the characteristics of the evolution equations satisfied by the 
constraint equations. Thus it is possible to understand that the change in boundary conditions 
induced by the change on the direction of the characteristics is just the one needed to ensure the 
correct preservation of the constraint equations under evolution. Among all possible evolution 
system we take a specific one which is both the simplest in terms of computations needed and 
also on dealing with constraint preservation. We then look at the residual gauge freedom of 
the Eddington-Finkelstein coordinates, that freedom is the only dynamics which is present in 
the vacuum case, and is the mode that can generate instabilities. After that we discuss the 
initial-boundary value problem for the specific evolution system chosen and determine which 
fields must be given in order to have a unique evolution, in particular we discuss a case in which 
the characteristics change at the inner boundary according to whether that boundary is inside 
or outside the horizon. We then discuss the instabilities of the system, due to the simplicity of 
the evolution system chosen this task is trivial and we can display explicitly the instability. We 
can decide in which cases there will be problems and the type of them. In particular near the 
stationary regime, that is when most of the scalar field has dissipated away or fallen into the black 
hole the system can be made stable using simple initial-boundary values for the gauge fields. 

Having completed the analytic study of the equations in section ^ we discuss the numerical 
methods employed, look at the convergence tests and constraint preservation under evolution for 
several types of initial-boundary data sets. 

In the next section we comment in several results obtained by the simulations, they do not 
pretend to be novel or interesting, but are done just to confirm the good properties of the evolution 
system chosen. They include ringing and tail studies, the relation mass-flux for boundary data 
and the relation total mass black hole mass for different types of collapse. 

2 The equations 

We consider the Eddington-Finkelstein coordinate system (|^), the time dependent spherically 
symmetric metric is 

ds^ = (-a2 + a2 dt^ + 2a^p dt dr + dr^ + r^b^ dft^ 

with a, 5, a and f3 being functions of r and t, and dfl^ the metric of the unit-sphere. The areal 
locking condition {b — 0) and the condition that the vector dt — dr be null is achieved choosing 
the lapse and the shift as 

_ raK^e _ a 

^ ~ l + raK% l + raK^e' 

where g is the corresponding coordinate component of the extrinsic curvature of the constant 
time hypersurfaces. Moreover this choice corresponds to the choice 6=1, and with it the metric 
becomes 

ds^ = a^i{2(3 -l)dt + dr) {dt + dr) + r^ dfl'^. 
The full set of evolutions and constraint equations for the geometric variables (a, K^e, K^r) are ^ 

dta = dr{a(3)-a^(l-l3)K\ 

^ Notice some corrections with respect to the equations in also notice the difference on the choice of 
gravitational constant. 
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respectively. While the wave equation for the massless scalar field becomes the first-order (in 
time) system 



(3) 



dtW = (r2(/3n + (i -/?)$) 



where n = 9t0, and <I> = dr(t>- 

In order to obtain first order evolution equations, we perform the change of variables (/, g, h, x+, X-) 

f = K\ar + l g=j h = {^drf + drU - a^K^). 



X+ = 



(n + $)r 



X- = 



(n- $)r 



/-I 



with this change the shift becomes f3 = — j — and 

ds^ =g{{f-2)dt + f dr) {dr + dt) + dl^^. 
Then the full set of evolution and constraint equations for the metric variables (/, g, h) are 



dtf 
dt9 
dth 



df ^ 9-2 + f 
dr r 
g8f^f-2d^^^^^_g_ig-2 + f) 



f dr f dr 



(4) 



dh 1 df a-2dg {g-f)h {g-2 + f)^^.x++X 



dr 2fr dr 2 fhr dr fr 



2/r2 



+ 2(- 



fr 



and 



Cf 



^-2if-2)h+l±l^-S = 
dr r jr 



dg 2g{hfr-xl+X-) 



~ dr fr 
while the wave equation for the massless scalar field becomes 
dx+ if - 2)x- 







(5) 



dtX+ = 
dtX- = 



dr fr 
f-2dx- , 2x 



f dr f^r 



+ ^i2r{f~2)h + 2-g-f) + 



8{x-)xl X+ 



^ (6) 
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where we have aheady used one of the constraint equations. 

Clearly the evolution equations for the scalar field are symmetric hyperbolic as a subsystem 
and their characteristics are fixed, they correspond to the two null directions. On the other hand 
the set of evolution equations for the metric coefficients is not unique, for we can add to them 
terms proportional to the constraint equations and change some of their characteristics at will. To 
study this possibilities we add to the evolution equations for /, g, and h respectively the following 
terms {Kfj - l)Cf + Kf.Cg, {K,f + j)Cf + {K,, - fj^)C,, {K^f - §^)Cf + [Kn, + ^^)C,. 
We get a system with the following principal part matrix, |^ 

Kff Kfg 

Kgf Kgg 

Khf Khg 1 

The characteristics of this system are given by the one forms (1, 1), (1, A+), and (1, A_), where 
A± = (if// + ± ^{Kff~KggY+AKfgKgf)/2. 

The first one corresponds to the propagation of h and it does so along the null incoming direction. 
It is the only fix characteristic. The other two can take any values, for instance if we take 
Kfg — Kgf — then the characteristics are (l,if//), and {l,Kgg), thus choosing these two 
values we can prescribe any propagation direction we please, they can be incoming, or outgoing 
time-like, null or even space-like. Even more, choosing Kff — Kgg — 0, Kfg = —Kgf ^ the 
characteristics become imaginary and so the system is not even hyperbolic! So we see that the 
different ways of writing Einstein's equations are only equivalent in the sense that solutions of 
one are solutions of the others, but some of their evolution equations do not even have a well 
posed initial value problem. We have tried to evolve numerically a system as the last one (with 
Kfg = —Kgf = 1) and as expected it is completely unstable, giving an explosive blow up of higher 
wave number modes. To avoid the fact that no good boundary conditions are known for non- 
symmetric hyperbolic systems these equations were tested using periodic boundary conditions 
and the initial data did not satisfy the constraint equations. But the instability does not seems to 
be sensitive to the failure of the constraint equations to be fulfilled nor to the particular boundary 
values chosen for those instabilities appear locally at every point of the evolution hypersurface. 

Thus, to start with, we must choose coefficients such that the system is hyperbolic. Since 
the coefficients Khf^ Khg do not alter the characteristics we shall take them to vanish, thus 
simplifying a bit the numerical computations. The sign of the characteristics at the boundary 
points determine whether or not values for the different fields must be given or not as boundary 
conditions. But not all these values are free, for they must be consistent with constraint propa- 
gation. In order for the constraints to remain enforced along evolution we must make sure that 
the resulting boundary conditions for the evolution equations obeyed for the constraints ensure 
that they have a unique solution, and it is the trivial one. It follows trivially from the form of 
the constraints that the characteristic matrix for the constraints evolution equations is: 

( Kff Kfg \ 

\ Kgf Kgg ) 

So that the propagation directions for them coincide with the propagation directions of the 
pair {f,g). This fact seems to imply that the freedom in choosing part of the pair {f,g) at 
the boundaries is just the one needed for getting maximally dissipative boundary conditions 
for the constraint system, thus making sure they propagate correctly. Unfortunately even for 
this simple case a rigorous proof of this has proven be difficult ^, and the resulting conditions 

^Recall that the principal part matrix is the one form by the coefficients of the derivative terms appearing on 
the equations. 

^The main difficulty arises when the K coefficients depend on the metric components 
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are complicated to implement numerically. Thus it seems that the best one can do retaining 
simplicity is to choose coefficients so that the characteristics for the pair (/, g) are incoming or 
tangent in both boundaries. This way no boundary condition is needed for the pair (/, g) and 
therefore no boundary condition must be satisfied by the constraints, thus, there is no danger 
that the constraint would cease to be satisfied during evolution. The simplest case of all is when 
both are tangent, and this is the case for instance if we take all coefficients to vanish. With this 
choice the evolution system reduces to. 



dtf - 2(/-2)/i 



4x^ 



rf 

2gxl _^ 2ff(/ - 2)xl ,_, 
dtg = 2gh-^ + — (7) 



dth 



And similarly the evolution equations for the constraints are just ordinary differential equa- 
tions, 

dtCf = ^{h fr -2x1) Cf (8) 



= ^{fxl-fxl+xl)Cg + ^{hrf-fxl-2xl + fX-)Cg (9) 



Since they are homogeneous ordinary differential equations the constraints will remain to be 
satisfied along evolution if they are initially satisfied. No extra condition is need at (nor can be 
imposed on) the boundaries. This is the system we shall study in detail from now on. Note that 
this system not only is symmetric hyperbolic, but it is also diagonal. In fact we reach such a 
simple system when looking for variables which would give diagonal systems. In figure ^ we show 
all characteristics of the chosen system. 
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2.1 The Gauge Freedom 



The chosen lapse and shift functions depend on the dynamical variables and therefore do not fix 
completely the gauge. To display this freedom we introduce the null coordinate -u = + r) ^. In 
the (u, r) coordinates the metric becomes 



For Schwarzschild the choice g — 1, f — 2 = —(1 — 2M/r) gives the standard form of the metric. 
If we choose some other coordinate u — u{u) the functional form remains the same if g = g/^ 



The invariant ratio is related to the mass function, indeed, 



is the usual form of the local mass for spherical symmetric space-times. 

Another invariant is the surface where / takes the value 2, this can be seen to be an apparent 
horizon. 

From the first constraint equation, (|^), it follows that under this gauge transformation, in the 
vacuum case. 



To fix completely the gauge we must uniquely define the coordinate u, this can be done by 
choosing their values at the initial time slice and at the outer boundary of the integration region. 
We shall do this by choosing the values for h at those points. 

2.2 Initial Boundary Value Problem 

We are interested in performing a time evolution between two fixed (areal) radiae, Rm, Rout, 
starting at a given initial surface, t = Iq. Thus we have to see what initial-boundary values can 
be prescribed. 

With the standard choice of equations (||,|^) it is clear that we need to give initial values for 
the fields (/, g, h, x+, X-)- Since they must satisfy the constraint not all initial data can be freely 
specified. We have taken a set of initial data sets for which the constraints are easily solved. We 
give an arbitrary initial value for x+ and solve 



ds^ = 9i{f - '2)du + 2dr)du + dVt^ , 



and f^^l^. 



du 




or 



dm 



Then the following choice is a solution of both 



constraint equations, 



/ 
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0, 
1, 

1 + 2m/r 



(10) 



h 




(11) 



*We thanks M. Tiglio for pointing to us this way to proceed. 
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In our numerical simulations we have taken x+ to be a Gaussian or a power of r times an 
exponential decay, so that the mass integration could be done exactly. The integration constant 
for m, the mass at can take any value, in general we have taken either |^ ™|Ki„ = 1 or 

We now discuss the boundary conditions (see figure (^. If the inner boundary is inside the 
apparent horizon (which is located at / = 2) there is no incoming characteristic at it and so no 
boundary condition is needed (nor can be prescribed) there. In fact in this case this boundary 
is actually space-like. If the inner boundary is outside the horizon there is one characteristic 
incoming into the integration region, namely the one for x_, so we have to prescribe some value 
for it. To allow for an apparent horizon marching across the boundary the numerical code checks 
for the value of / at the boundary and, depending whether its value is smaller or bigger than 
2, imposes a boundary condition for x-- When needed we have normally used a null incoming 
radiation condition, {x-\Rin — 0) but it is clear than in that case there is no natural boundary 
condition which can mimic the physical collapse one is trying to simulate. 

At the outer boundary, and provided the horizon is towards the inside of that boundary 
(/Iflouter < 2), there are two incoming characteristics, one corresponding to x+ the other to h. 
We can prescribe arbitrary values for them there, and we have performed several runs prescribing 
either incoming scalar field radiation or pure gauge modes. 

With these initial-boundary value conditions the problem is well posed and so, given smooth 
data we obtain a smooth solution valid for a finite time interval. 



2.3 Analytic Instabilities 

The evolution system is well posed, so there are no instabilities of the explosive type, that is 
those whose rate grows unbounded with the frequency, rather the expected instabilities have an 
exponent with bounded positive real part. These instabilities essentially affect the pair (/, g) 
which controls the characteristics of the system, if they are present, and act for a long enough 
time they can ruin the hyperbolicity of the system -by making some of the propagation speeds 
to diverge- and so its well posedness. But much before that happens one of the characteristic 
speeds grows so much that the system becomes numerically unstable due to violations of the 
Courant-Friedrich stability condition. 

For simplicity we analyze the instability for the vacuum case. In that case the equations are: 



dtf = 2{f-2)h 

dtg = 2gh (12) 
or 

Given initial data for (/, g, h) and boundary data at the outer boundary for /i, we can first 
solve for h, h{t, r) = hi,{t + r — Rout) and then for (/, g). We have. 



f{t,r)-2 = (/o(r)-2)e'/o'*''(*+'-^°"*^''* 
g(t,r) = go(r)e Jo 



Thus we see that g remains always positive, and the horizon -the radius at which f{t, r) — 2 
remains fixed. But the equal time surfaces tilt, and so if h takes big positive values, /, and g grow 
exponentially and we have an analytical instability, which in numerical simulations would result 
in a numerical one if no appropriate methods are used. If a scalar field is present, no matter how 
small, the tilt of the equal time surfaces would increase the propagation speed of the outgoing 

^ Due to the scaling properties of the system there is no loss of generality in choosing only just these two values. 
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modes, v — '^-j^- At points inside the horizon (/o — 2 > 0) and so there / — 2 grows exponentially. 
The propagation speed v goes to 1. At points outside the horizon (/o — 2 < 0) and so there / 
eventually becomes negative. Notice that when / vanish the equal time hypersurface becomes 
null and so the time evolution of all outgoing modes (in this case only the one from the scalar 
field) is ill posed in this gauge. 

On the other hand if h becomes large but negative, / ^ 2 and now the hypersurface of 
constant r approach null surfaces. In that case the evolution still can be done, but the loss of 
accuracy is very important, even for moderate values of h. 

We need therefore to keep h moderately small during evolution and propagating inwards 
without any important residual part due to the numerical method. This can be achieved in many 
cases by just prescribing the initial-boundary conditions for h so that it starts small. There are 
nevertheless situations for which h becomes big and the simulation runs into the instabilities, 
this is so, for instance, if one takes boundary data for the incoming component of the scalar field 
with a very big amplitude (final masses of order 10). In that case the h generated along the 
evolution of the scalar field is important and the systems becomes unstable. On the other hand, 
for situations describing the final stages of collapse, namely after most of the scalar field has fallen 
into the black hole and the geometry has settled in a near stationary regime the instabilities do 
not show up. We have tested that with evolutions lasting for several thousand masses. 

3 Numerical Methods 

To perform the simulations we use a standard fourth order method of lines. The space derivatives 
discretization is done with standard fourth order (we have also used second order discretizations 
with similar results) discretized derivatives which are one sided at the boundaries. We have used 
symmetric and non-symmetric, but full fourth order, ones. The symmetric ones, pO| , |ll|, are 
guaranteed to be stable, but are only second or third order accurate at the last points of the 
grid near the boundary. This accounts, for instance, for spikes at the error on the constraint 
propagation check, some of these errors sometimes propagate inwards, but stay bounded and 
diminish as the grid spacing is decreased. The non-symmetric discretizations used were full 
fourth order (but with bigger errors near the boundary) and they give smaller errors on constraints 
evolution near the boundaries. But there is no guaranty that the scheme would be stable. Most 
of the calculations, including all the results and test given in this paper, were done with the 
symmetric ones. Their symmetry implies they are norm preserving and so they do not have any 
dissipation. Fourth order centered differences approximations introduce hight frequency modes 
with backwards propagations speeds (group velocities) which are 5/3 bigger than the exact speeds. 
These modes become important when analyzing the tail behavior of the scalar field, for there the 
values of the scalar field become so small that the contribution from them becomes important. 
The higher than light propagation of them implies that boundary effects take place before they 
should. These high frequency modes become smaller as the grid size is decreased, thus tail decay 
studies are possible if enough grid points are included, but this should be very costly in higher 
dimensions. For these studies it seems better to use centered difference schemes of second order, 
for they have a maximum group velocity of one. 

Time discretization is done using the standard fourth order Runge-Kutta method. 

We have done three types of simulations: The first type of simulation has no-incoming bound- 
ary conditions {x+ = /i = at the outer boundary) and an initial data satisfying (p^. One 
starts with a solution representing an inner black hole (usually of unit mass) and certain amount 
of inf ailing scalar field. The evolution should result in a bigger black hole and some amount of 
scalar field going outwards. Thus a ringing and latter a power law tail should be observed in the 
outgoing component of the scalar field, x-- 

The second type are vacuum solutions representing pure gauge dynamics. They have Schwarz- 
schild or Minkowski initial data and the dynamics is created by allowing incoming gauge {h\u^^^ = 
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h^{t)) to enter the outer boundary. The evolution should have some dynamics and latter on, if 
the incoming gauge has a finite duration, should relax into a static black hole in a different gauge, 
through the simulation the mass should stay constant, as should the horizon position. 

The third type of data starts with initial data corresponding to a Schwarzschild black hole of 
unit mass and latter on the evolution an incoming scalar field mode is injected into the system 
from the outer boundary during a finite amount of time (x+ 1 flout = X+(^)- The evolution should 
create a bigger black hole and some outgoing radiation should be seen. 

We shall discuss the results for these three type of simulations in the next section. 

3.1 Convergence Test 

We have perform convergence tests for the three types of simulations described above. In general 
it is difficult to analyst the convergence order in detail due to the fact that the space discrete 
derivative we use is not of the same order near the boundary (third or second order) rather than 
fourth. This gives rise to Q values different than the corresponding to fourth order (Q=16). 
Typically after boundary effects become important (that is when the fields reach the boundary) 
we get values over 12, but some times they drop to about 4 when the bulk of the solution is near 
the boundary. In figure: ||, and || we give values for the Q factors obtained using the and 
L°° norms respectively, for a run with homogeneous boundary conditions and a Gaussian peak 
of the ingoing component of the scalar field midway in the integration region. This run has an 
initial black hole mass Mjniuai = 1, a total mass, Mxatai = 8.90 and a final black hole mass 
of Mpinai — 8.67 so it is well into the nonlinear regime. The peak observed at the end of the 
integration period is due to boundary effects, after that the values of the different Q's go to values 
above 12, but never return to 16. 




3.2 Constraint Evolution 

The constraint equations are not solved for in any step of our simulations, and although they 
are satisfied along evolution for the exact equations here we have to check this happens at the 
discrete level. The figures (||-§1) show the L°° and norm of the constraint expressions along 
evolution of different type of data set used. For type I data we had: MjniUai = 1, Mxotai = 1-15, 
Mpinai = 1.12. For the type II data we introduced at the outer boundary 10 cycles of gauge 
mode of frequency 27r, starting at t — 10. The black hole mass was of course constant and had 
the value 1. For the type III data we introduced at the outer boundary 10 cycles of the x+ mode 
with frequency 27r, starting at t = 10. We had: Mjnitiai = 1, Mxotai = 2.6, Mpinai = 2.4. 
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Figure 4: Log of the L°° norm of C/ of type 
I data 
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Figure 6: Log of the norm of C/ of type 
II data 
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Figure 8: Log of the norm of Cf of type 
III data 
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Figure 5: Log of the norm of Cg of type 
I data 
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Figure 7: Log of the norm of Cg of type 
II data 
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Figure 9: Log of the norm of Cg of type 
III data 
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It is clear that the constraint quantities remain bounded along evolution and that they go 
to zero as the grid size diminishes. Most of the contribution to the L°° norm comes from the 
boundary, where the derivative operator used is only second or third order. These peaks near de 
boundary can be reduced to the same level than the rest if one uses a discrete derivative operators 
which are fourth order at all grid points. 

If one looks at the ratios C*//^, ^g/^- which in some sense measure better the failure of the 
constraint to be satisfied one sees that those ratios do not change so much along evolution. For 
Cf the change is of only about 7%, and are of the order of lO^"*. For Cg, since g is chosen initially 
to be zero, the change at the beginning is rather big, but after a transient it gets to a plateau of 
the order of 10"'*. 



4 Results 

4.1 Decay of the scalar field, ringing and power law decay 

We reproduce standard results on ringing and tail (power law decay) for the first type of simu- 
lation. In figure ( |l0|) it is plotted the log of the absolute value of X- at r = 25.5 with respect 
to time. It is a run where the outer boundary is at Router — 6000, the space grid spacing was 
h/M = 0.0234, where we have used the initial mass to make the ratio, footnotein most of the 
runs k/h k, 2, this does not violate the Courant-Friedrichs-Lewi condition because we use Runge- 
Kutta for time evolution. The initial scalar field is given by a Gaussian x+ pulse centered at 
r = 10, giving: Initial black hole mass Mjnitiai = 1, total mass MTotai — 2.56, final black hole 
mass Mpinai = 2.45. At the beginning the characteristic ringing is seen and then a power law 

d Inf V ) 

decay of the type (this can be seen best plotting the local decay power t — as defined in 

1^, |l^, figure (^l])). The power law decay computed t~^-°^ (at t — 7200), agrees quite well from 
the one expected for linear perturbation theory, namely See |l^. In the plot we also 
show a lower resolution {h/M = 0.0468) computation where it is clear that an instability starts 
to appear and become important at about half the time. This instability diminishes when the 
grid resolution is incremented, but still seems to be present at longer times. At the resolution 
used a feature departing form the power law decay is seeing after t = 7200, this corresponds to 
an higher frequency mode propagating backwards at a speed 1.66 the speed of light. These high 
frequency modes are common on all centered fourth order methods. Their amplitude decreases 
with increased precision. The time of arrival of that feature implies it originates near the inner 
boundary probably due to imperfect ingoing boundary conditions. After that feature, we again 
see (not shown in the figure) a even bigger feature, this one propagates at the correct speed of 
the problem, but its presence at that time implies that it is generated by the one traveling at 
higher speed when it enters the region near the outer boundary. 

The possibility of giving boundary conditions which would automatically satisfy the constraint 
equations allow as to perform several interesting runs. In particular we made a couple of runs of 
type III with 

= A+{1 - cos(w+t)) t G [10, 20] 

with j4_|- = l/27r, and in one case lo — 27r, in the other uj = 27r/10. For these two boundary data 
we have, 

/.20 

/ xldt = IbAl 

J 10 

Since the energy fiux is given by 

dt g^^+ P 
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Figure 10: Decay of X- at r = 25.5 Figure 11: Power law decay of X- at r = 25.5 

and ^ « If x+, It ~ J?x\, %^ ~ ^X+- We see that to a 1% accuracy the flux of both 
solutions should coincide and so the black hole masses should be similar. Even more, the ringing 
of both solutions should be very similar, for it mostly depends on the geometry. In figure (|l^) we 
see the mass of one of the simulations and the mass difference between both simulations multiplied 
by 5. We see that most of the difference is bounded to the region where the wave of scalar field 
is moving, leaving behind the same geometry. In figure (|l^) we see one of the ringings (the value 
of X- at r = 50) and the difference between the ringing of both solutions augmented by a factor 
10. Again we see that most of the difference is at the moment where the two different wave fronts 
pass the point and then on the precise location of the maxima of the ringing. 

-1 ^ , ■ ■% 




Figure 12: Mass of type III run and mass difference between runs augmented 5 times. 



4.2 Relation between the initial mass and the black hole mass 

We have tested the relation between the mass gap given by the total initial mass minus the final 
black hole mass as a function of the total initial mass for type I and type III initial-boundary 
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Figure 13: Values of X- at r = 50 for a type III run, bellow the difference between runs augmented 
10 times. 



data sets. In figures (14, y_5|) we display the results. In both cases we start with a black hole of 
mass unit, and put our inner boundary inside its horizon. 
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Figure 14: AM vs M for type I data, AM = 

0.21 * (M - 1)0-86 g-0.55(Af-l)i « 



Figure 15: AM vs M for type III data 



5 Conclusions 



It is difficult to see whether some of the ideas presented in this work can be extended to the full 
three dimensional case, where many more variables are present and where there is no preferred 
center in which to anchor a gauge as the one used here. But the model is so simple that perhaps 
the observations made in solving it can shed some light into the more difficult problem. Among the 
observations we have the following: If the gauge prescription docs not fix completely the gauge 
it is expected there would be gauge modes propagating. In this case there was just one gauge 
freedom left (the value of h at the initial surface and at the outer boundary), nevertheless there 
were three equations which acquired a non trivial propagation, h) Gauge modes, by being non- 
physical, can have singular behavior although the geometry given by the solution can be regular. 
It is clear in the model studied that the different components of the metric grow exponentially 
while the solution is still Schwarzschild. Thus, it seems to be necessary to isolate and keep under 
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control the potentially unstable gauge modes. In our case this was done by choosing variables for 
which it was clear where the instability was, and so, choosing conveniently the initial-boundary 
data we could keep it small for many cases, although not always, c) The choice of gauge modes 
can also complicate the boundary value problem. Ideally some of the fields should acquire values 
at the boundary so that the constraint propagate correctly. In the model it can be seen that if 
part of the data for the pair (/, g) must be given at the boundary (according to the values chosen 
for the coefficients {Kff,Kfg,Kgg,Kgf)) that part must satisfy some evolution equation along 
the boundary, so actually there is no freedom left. But other gauge quantities, like h, can take any 
value. Thus, if the variables are not chosen appropriately it could become very cumbersome to 
find the correct boundary conditions which would, at the same time keep the constraint equations 
being satisfied and the gauge instabilities under control. Note that the procedure we used, getting 
some equations which are intrinsic to the boundary for some of the variables, is similar to the 
one used in JTsf , the only case where well posedness of the initial-boundary value problem in full 
general relativity has been asserted. 
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